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We derive various sampling functions for multimode homodyne tomography with a single local 
oscillator. These functions allow us to sample multimode s-parametrized quasidistributions, density 
matrix elements in Fock basis, and s-ordered moments of arbitrary order directly from the measured 
quadrature statistics. The inevitable experimental losses can be compensated by proper modification 
of the sampling functions. Results of Monte Carlo simulations for squeezed three-mode state are 
reported and the feasibility of reconstruction of the three-mode Q-function and s-ordered moments 
from 10 7 sampled data is demonstrated. 
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I. INTRODUCTION 

Recent development of quantum-state reconstruction 
methods has made it possible to completely reconstruct 
an unknown state of a quantum mechanical system pro- 
vided that many identical copies of the state are avail- 
able. The method was pioneered in quantum optics, 
where an optical homodyne tomography was devised to 
reconstruct a quantum state of traveling electromagnetic 
field Other proposed techniques involved unbal- 

anced homodyning Q and cavity-field measurements by 
atomic probes Q. Quantum-state reconstruction proce- 
dures were also successfully applied to molecular vibra- 
tional state and motional quantum state of trapped 
ion || . 

Single-mode optical homodyne tomography is now a 
well-established technique. Based on balanced homodyne 
detection, the method seeks to reconstruct the quantum 
state from the statistics of the quadrature components 
of the signal mode. The standard experimental setup 
involves balanced lossless beam splitter, where the sig- 
nal is mixed with a strong coherent local oscillator (LO). 
Two photodetectors are placed at the beam splitter out- 
puts and the two photocurrents are subtracted, thereby 
removing LO fluctuations from the resulting signal. 

Wigner function can be obtained from the measured 
quadrature statistics by means of inverse Radon trans- 
form HJ. Once the Wi gner function is known, expecta- 
tion value of any operator can be be evaluated by averag- 
ing corresponding phase-space function over the Wigner 
quasidistribution. This strategy, however, is not optimal, 
because the experimental errors are amplified during nu- 
merical data processing and the final error can be very 
large. Fortunately, the detour via Wigner function can 
be avoided and density matrix elements in Fock basis can 
be directly reconstructed by averaging appropriate sam- 
plin g functions over the measured quadrature statistics 
J10|— |l5|| . The functions for sampling s-ordered moments 
were found in pq-p^|, and those allowing direct recon- 
struction of the exponential moments of quantum phase 



distributions were obtained in JljJ^Cj (for a review, see 
|^l| , p2| ). The problems with inverse Radon transform can 
be avoided by reconstructing smoothed Wigner functions 
p3| . The sampling is a simple and straightforward lin- 
ear operation which can be in principle performed in real 
time during experiment. We note that, besides linear 
sampling procedures, reconstruction strategies based on 
maximum likelihood estimation and maximum en- 
tropy principle |25j ] have been proposed. 

Recently, increasing attention has been devoted to 
multimode homodyne tomography |prj|-|32|| , because some 
of the most interesting quantum mechanical phenomena 
stem from correlations between several degrees of free- 
dom. Let us mention just the EPR paradox and vio- 
lation of Bell's inequalities |||]. Quantized electromag- 
netic field is one of the most suitable systems for thor- 
ough investigation and exploitation of these phenomena. 
For example, entangled signal and idler photons can be 
routinely prepared by means of spontaneous paramet- 
ric down-conversion p3]. The entangled photon pairs 
play crucial role in certain quantum state teleportation 
schemes [|35| and quantum cryptography setups |]36f . 

Multimode extension of optical homodyne tomography 
is straightforward. One can introduce a separate LO and 
homodyne detector for each mode of interest and mea- 
sure joint multimode quadrature distribution. In this 
case the sampling functions developed for single-mode to- 
mography can be immediately employed. Very recently, 
this approach has been used to measure the joint photon 
number statistics of two-mode squeezed state prepared 
in a nondegenerate optical parametric amplifier 

However, the requirement of specific homodyne detec- 
tor for each mode complicates the experiment. It would 
often be much more feasible to use only one homodyne 
detector and one LO. In such experiment, a distribution 
of one quadrature X, which is a linear superposition of 
N single-mode quadratures, is measured. The knowledge 
of the probability distribution of all distinct quadratures 
X provides a complete information on the multimode 
quantum state. In particular, two-mode tomography 
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with single homodyne detector was discussed extensively. 
The functions for sampling density matrix elements in 
Fock basis were found in |27],^,|3l|] , and those allowing 
direct reconstruction of two-mode correlation functions 
were obtained in p9|-|3l| . Recently, a general multimode 
homodyne tomography with a single LO was considered 
and the sampling functions for density matrix elements 
were expressed in terms of integrals |32j]. In this paper 
we shall derive various important sampling functions for 
multimode homodyne tomography with a single LO. All 
functions are expressed in analytical form. Imperfect de- 
tection is considered and it is shown that the losses can 
be compensated by proper modification (rescaling) of the 
sampling functions. 
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FIG. 1. Measurement of internal quantum correlations of 
optical pulses [29]. The signal pulse and a train of strong 
local-oscillator (LO) pulses that are short compared to the 
signal pulse are mixed at a 50%:50% beam splitter (BS) and 
the two photocurrents measured by photodetectors (PD) are 
subtracted. The train of LO pulses is prepared interferomet- 
rically, thereby allowing one to control the pulse distances, 
relative phases and intensities, as symbolically denoted by 
U(O it il>i). 
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FIG. 2. Optical homodyne tomography of sin- 

gle-frequency multimode optical field. The signal modes 
a\ , . ■ ■ , ajv feed the input of N-port interferometr which pre- 
pares the mode A at one of its outputs. Subsequently, the 
distribution of quadrature X is measured by means of stan- 
dard single-mode homodyne detection. 



The paper is organized as follows. In Sec. II we ad- 
dress the reconstruction of multimode smoothed Wigner 
functions. The results are then applied in Sec. Ill to 
find the sampling functions for density matrix elements 
in Fock basis. The reconstruction of s-ordered moments 
of the field operators is discussed in Sec. IV. The results 
of Monte Carlo simulations of multimode homodyne to- 
mography are reported in Sec. V. Finally, Section VI 
contains conclusions. 



II. SAMPLING FUNCTIONS FOR 
^-PARAMETRIZED QUASIDISTRIBUTIONS 

In multimode homodyne tomography with a single lo- 
cal oscillator one measures a probability distribution of 
the quadrature 



X.-tfA + At), 



(i) 



where the operator A is a linear superposition of annihi- 
lation operators a,j of N signal modes, 



A'* 



'3 3 ■ 



(2) 



The complex coefficients z/ fulfill normalization condition 



N 

Ei 

3=1 



(3) 



which ensures validity of standard commutation relation 
[A^l = 1 for the operator A. Two examples of ex- 
perimental setups, where the statistics of quadrature X 
are measured, are given in Figs. 1 and 2. Multimode 
homodyne tomography can be employed to investigate 
ultrafast internal quantum correlations of optical pulses 
p9| , p0[ , see Fig. 1. A train of N strong LO pulses is 
used to select a set of N nonmonochromatic modes from 
the signal pulse. The modes a,j are determined by posi- 
tions and shapes of the LO pulses and the correlations 
of the signal pulse are probed in terms of these modes. 
Figure 2 illustrates a scheme for homodyne tomography 
of single-frequency multimode field. The desired super- 
position A is prepared in Alport interferometr and then 
it enters homodyne detector. This setup can be used 
e.g. for measurement of a polarization state of an optical 
field [j37|| . The modes a\ and ai then correspond to two 
orthogonal linear polarizations. The two-mode unitary 
transformations U(8, ip) leading to superpositions (|J) can 
be performed with the help of two phase shifters and a 
polarizing beam splitter |]37| . A common feature of the 
experimental setups shown in Figs. 1 and 2 is that only 
one balanced homodyne detector is needed. 



If the statistics w(X; {zj 
known for all {zj} fulfilling 



) of the quadrature X are 
t) , then we have a complete 
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knowledge of the quantum state of the multimode light 
field and all quantities of interest, such as various qua- 
sidistributions, density matrix elements, and s-ordered 
moments, can be unambiguously determined from the 
distributions w(X; {zj}). 



A. Sampling of the smoothed Wigner functions 

Let us begin with reconstruction of the multimode 
s-parametrized quasidistributions. It is convenient to 
work in the hyperspherical coordinates. The points {zj} 
lie on a surface of 2A-dimensional unit sphere and we 
parametrize them as |32] 



where 



3-1 



Uj (6) = cos 9j Y[ sin t , j < N 



AT-l 

un(0) = Y[ sin ^' 
1=1 



(4) 

(5) 
(6) 



and 



^•e[0,27r], j = l,...,N, 
^•G[0,7r/2], j = l,...,N-l. 

To simplify the notation, we define 9 — . . . ,0jy_i) 
and ij> = (ipx, . . .,ipjf)- 

Multimode characteristic function corresponding to s- 
ordering of the field operators is defined as |58| , 

C.dPj}) = (j[ exp Q S |/3,| 2 + ja ] ~ /3*a^j J , (7) 

where () denotes quantum mechanical average. Let us 
compare the exponent on the right-hand side of Eq. (|t]) 
with the quadrature X({zj}) = X(0,ip). We can see 
that C s ({Pj}) is proportional to characteristic function 
of the quadrature distribution, 

/>oo 

C s ({&}) = e^ 2 / 2 / dXe^ rX w(X;e,i>), (8) 



where (3j = iruj(0) exp(?/0j) and r > is radial variable, 



N 



Multimode quasidistribution W s ({%'}) is a Fourier 
transform of the characteristic function C s ({(3j}), 

N 

W.({<*j}) = -5F / Cs({Pj}) II d2 & ePfr-W. (9) 

J 3 = 1 

We rewrite this integral in the hyperspherical coordi- 
nates. We shall integrate over the angles 6j, phases ipj, 
and radius r. It is convenient to introduce dfl, 



JV-l N 

dn = g(9) J[ d0ifj#;> 

1=1 ]=1 

where the prefactor 

JV-l 

g{6)= cos^CsineO 2 ^- - 1 
i=i 



(10) 



(11) 



stems from the Jacobian of coordinate transformation. 
We substitute the characteristic function (||) into (||) and 
after some algebra we arrive at 

-i />oo /> poo 

W s {{a 3 }) = — / dr dn dXr™- 1 
~ Jo Jn J-oo 

where we have introduced a c-number quadrature 
1 N 

X({ aj },d,^) = (uj^e-^aj +c.c.) . (13) 

.7=1 



After changing the order of integration in (|12|), we find 
that 

W s ({a j }) = [ dn f dXw(X;0,ip) 

Jn J-oo 

xS N (x-X({ aj },9,iJ>);s) , (14) 
where the sampling function Sn reads 

S n& s ) = ^nJ q dre^/V^r 3 *- 1 , s<0. 

(15) 

This expression can be further simplified. The quasidis- 
tributions Ws({aj}) as well as the quadrature distribu- 
tions u>(A; 0, ip) are real functions. The sampling func- 
tion Sn has to be real and only real part of the above 
integral should be considered. The imaginary part of Sn 
is a null function whose average over any physical quadra- 
ture distribution w(X; 0, tp) is zero. Thus we can replace 
exp(zv / 2^) by cos(-\/2V£) in Eq. (jl^)- The integration 
over r can be easily carried out and yields a confluent 
hypergeometric function, 

2 N ~ 1 (N-1V ( 1 £ 2 \ 



The function Sat depends on ay, A, 0, and i/> only 
through a specific combination £ = A — A. 

The parameter s must be negative because the inte- 
gral (|l5|) would diverge otherwise. This implies that 
only smoothed Wigner functions corresponding to s < 
can be directly sampled from homodyne statistics. The 
confluent hypergeometric functions can be expressed in 
terms of the error function of the imaginary argument 
erfi(x). It holds that 
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FIG. 3. Sampling functions Sn(£, — 1) for the Husimi Q-function of iV-mode optical field. 



$ ( 1, -; -x 2 ) = 1 - y/nxe ^erfi(x), 



, - i , A -'- A 



(-1) 



which allows for an easy determination of the required 
sampling function. 

Our results form a multimode generalization of the 
single-mode relations obtained by Vogel and Riskcn |{| 
and by Richter p3| . Notice also, that D'Ariano et 
al. gave explicit formula for sampling function of two- 
mode Husimi quasidistribution p^| . Several functions 
5jv(£;— 1) are plotted in Fig. 3. The number of oscil- 
lations of 5at(£;s) increases with increasing N and the 
sampling functions are bounded, Sn — > as |£| — * oo. 



B. Imperfect detection and loss-compensating 
sampling functions 

The sampling functions ( |l6| ) would yield correct re- 
sults only in the ideal case of unit detection efficiency. 
In a realistic experiment, losses are inevitable and the 
overall detection efficiency rj is lower than 1. The losses 
can be modeled as a mixing of the signal mode with a 
vacuum on a beam splitter. The detected quadrature X' 
is thus a superposition of the original quadrature X and 
a vacuum-state quadrature X vac p9|, 



X' = ^X + v/T^Xv 



(17) 



With the help of (fL7|) one can find a simple relation be- 
tween the characteristic functions of the quadratures X 
and X', 



exp (i\/2rX^ ^ = exp ^ - V 2 



cxp j i \j —rX' 



(18) 



Inserting formula ([Ts| ) into (H) and repeating the steps 
leading to Eq. (|l6|) one finds that the replacements 



X 



X 



(19) 



are necessary and sufficient in Eq. ( |l4| ) to account for 
detection losses, 



S N (x,X;s >V j =S N (^-X;s ' ' ' 



V 



(20) 



The losses impose a new limit on the ordering parameter 
s because the modified ordering parameter s + (1 — rj)/r) 
must be negative, 



^ 1 - '? _ 
s < = s v . 



(21) 



Only smoothed Wigner functions with s < s v can be 
reconstructed if losses are present in the experiment. 
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III. DENSITY MATRIX ELEMENTS 

In this section we briefly address the sampling of mul- 
timodc density matrix elements in the Fock state basis 
\{ni}) = \ni)\n 2 ) ■ ■ ■ \n N ), 

Pm* = {{m l }\p\{ni}) 1 (22) 

where m = mi, , . . , mjy and n = ni, . . . , njsr are vector 
indices used for notation simplicity. In tomography with 
single LO the matrix elements p mn can be reconstructed 
from the measured data according to 

p mn = [ dfl [ dXf mn (X,0,iP)w(X-9,4>). (23) 

The functions / mn were expressed in terms of integrals in 
Ref. [B2| . Well-known analytical formulas for single- mode 



sampling functions f mn involve products of regular and 
irregular eigenfunctions of the harmonic oscillator Hamil- 
tonian The two-mode functions / mim!; „ lB2 can 

be written as finite series of the confluent hypergeomet- 
ric functions pl|. Here we show how to derive analyti- 
cal expressions for arbitrary sampling functions / mn for 
generic TV-mode optical field. Our starting point shall be 
multimode Husimi Q-function, 

Q(R}) = 4v<R}HR}>> ( 24 ) 

where |{aj}) is multimode coherent state. When the den- 
sity operator p is expanded in Fock basis the Eq. (|24| ) 
takes the form 



3 3 r -\oci\ 



mi ,ni— rriN ,un—0 



i ^/m 3 \n 3 \ 



(25) 



From this expansion we can readily see that Husimi quasidistribution Q({aj}) = W-i({aj}) is a generating function 
of the density matrix elements in Fock basis, 



N i 

rr L 



N 



Q(K})II e|Q! 



;=i 



(26) 



It follows immediately that the sampling function for the Husimi quasidistribution is a generating function of the 
sampling functions / mn - This can be shown explicitly by inserting the expressions ( |l4| ) and (|2^) into Eq. (26) and 
comparing left- and right-hand sides of the resulting formula. We have 



N 



1 



in o _ rij 
3 



N 



S N [X, X({ aj }, 8, V); s = -1, rjj I] e l Q 'l 

1=1 



(27) 



a , —a* =0 



This expression is general, i.e. valid for any number of 
modes. The Q-function can be sampled only if the de- 
tection efficiency 77 > 0.5, c.f. Eq. (|2l|). This also limits 
the possibility of sampling the density matrix elements; 
the functions / mn exist only for 77 > 0.5. 

The dependence of / mn on phases ipj can be seen from 
Eq. (|27j ) even without going into explicit calculations. 
With the help of the substitution atj = jj exp{iipj) one 
obtains 



/v 



fmn{X,e,lP;f 1 )=F n 



3=1 



(28) 



moreover, F mn (X^6;ri) are real functions. Analytical 
formula for these so-called pattern functions -F mn can be 
derived if one inserts the sampling function Sn ( |l6| ) into 
( p7| ) and performs the necessary differentiations. After 
a tedious but straightforward calculation one finds that 



F mll can be written in terms of finite series of confluent 
hypergeometric functions, 



F mn {x,e-i 1 ) = 




where fij = max(m j -, rij), Vj — min(rrij, rij), 



5 



JV 



P»,u,k = ^2 Mj - v i + 2fc 3> 

3=1 



(30) 



and 



3jv(x, 2fc) = (-l) fe (7V + fc - 1)! $ (V + fc, i; -x 2 ^j . 
E N (x, 2k + 1) = 2x(-l) k (N + k)\ $ ( N + k + 1, |; -a; 2 



Notice an interesting analogy. The quantum state is 
uniquely and completely determined by its Husimi qua- 
sidistribution, which contains complete information on 
all density matrix elements p mn - Similarly, all the sam- 
pling functions for density matrix elements can be ob- 
tained from the sampling function Sn, which contains 
all information on / mn - 

A single-mode version of the formula ( ^6| ) was used in 
Jl0| to find the sampling functions for single-mode den- 
sity matrix elements. However, the sampling function 
Sn was not explicitly given in Jl0[ and the results were 
written in form of complicated series. Thus later different 
techniques have been adopted to calculate F mn [|l2] p"5| |. 
We emphasize that for N = 1 the Eq. ( p9f ) yields exactly 
the single-mode pattern functions F mn given in [ ^2|pT| 
and for TV = 2 we get the two-mode pattern functions 
derived in |2^] . The formula ( ^9|) is also suitable for inves- 
tigation of the asymptotic behavior. One simply inserts 
the asymptotic expansions of relevant confluent hyperge- 
ometric functions into ( p9| ) and extracts the asymptotic 
expansion of F mll . It turns out that all pattern functions 
are bounded and go to zero as \X\ — » oo. 



Finally we note that the functions F mn given by Q2S 
differ from those obtained in Ref. |32j because we have 
removed a superfluous imaginary part of Sn ■ Had we re- 
tained this imaginary part, we would have obtained the 
pattern functions derived in jj2). To see this, one can 
insert the integral representation ([l5]) into Eq. ( p7| ) and 
differentiate prior to the integration. One recovers the 
formula (13) of Ref. p§, 



F mn (X, 9- rj) = -4 fl \ ri H«i(»)] 

3=1 ' 



dre ^r* e i^rX r 2N-l 



xl[r»'-»'L%-»'lr 2 uf(8)}, (31) 



i=i 



where L"(x) denotes generalized Laguerre polynomial. 
A real part of the complex function (31) coincides with 
Eq. (g). 




FIG. 4. The functions F^(9) biorthogonal to G 2 k {6) in the 
interval [0, 7r]; m = solid line, m = 1 dashed line, m — 2 
dot-dashed line. 



IV. S-ORDERED MOMENTS 

A. Multimode sampling functions 

Here we consider direct sampling of the multimode s- 
ordered moments 



L/ mn — . . . u N u 1 ... u N / a . 



(32) 



We will follow an approach due to Opatrny et al. |29f | 
and generalize their results for two-mode homodyning 
to any number of modes. The quadrature distribution 
w(X; 6,iF) can be obtained from the joint distribution 
w(xi, . . . , xn] ip) of the single mode quadratures 



1 



according to 

w(x-e,tp) 



V2 (aje 



dxi 



(33) 



N 



:s\x-J2 Uj (e) 

3=1 




The moments (32) can be reconstructed from the joint 

s \ (m j +n j )/2 



quadrature statistics as follows 

N 



Cmn = J W(XI, ...,X N ]l/>)Y[ dx 3 dlpj (|) 

3=1 

fe) A ;;// ,. ;; , U ' " " . (35) 



rrij +n j 



where we integrate over N quadratures Xj <E (— oo, oo) 
and N phases ipj G [0, n]. H n (x) denotes customary Her- 
mite polynomial of variable x and 



K(m, n) 



m + n 



(36) 



G 



The multimode sampling function employed in Eq. ( pq ) 
is just a product of the appropriate single-mode sampling 
functions derived in 

We would like to link Cmn to the quadrature distribu- 
tion w(X; 6, ■0), 



cl s l = / dtl 



dXD mn (X;9,tP;s)w(X;6,^), (37) 



where 



JV-l 



1=1 



N 



d ^= n / Mi n / d ^ 



(38) 



Notice the definition interval of the phase variables, 
•ipj € [0, 7r]. The upper bound of integration over 9j is 
denoted by # max . The most straightforward choice would 
be, of course, to keep 9 niax — it/ 2 as in previous sections. 



As we shall see later, the choice # max = t can be more 
suitable. 

Following |^9| we shall look for the sampling function 
D m n in the factorized form, 



/ S \ M/2 / X 

L> mn (X,0,t/>;s) = M H M 



N-l 



K I II F m!+ni( l) 



N 

x J^J K(rrij,'i 

3 = 1 



1=1 



where M = X)ili( m ' + n i) an d ^m'+m (^/) are some yet 
undetermined functions. Inserting Eqs. (39) and (|34j ) 
into Eq. ( [37] ) and comparing the resulting expression 
with ( |3q ) we conclude that the following integral equa- 
tion must be fulfilled: 



(16 



N-l Hm 



1 N 

3=1 



N-l 



N 



n ^i+n^i) - W^+n, -4= 

1 = 1 3 = 1 



(40) 



We shall need the summation rule for Hermite polyno- 
mials, 



Hi (asi cos + z 2 sin 0) = ^ G z fc (0)£T fc (a:i) JTj_ fc (x 2 ), 



k=0 



where 



G\{6)= [ l k ){cos6) k {sin6) 1 -" 



(41) 



(42) 



If we use the summation rule (41) repeatedly we find that 



which is a multimode generalization of (|4l|). The prime 
denotes sum over all j%, . . . , j'jv meeting the constraint 

J2iLiJl = M i and 



N 



(44) 



The expansion ([13]) is inserted into Eq. (^) where the 
integration over 0/ should select the right sequence of the 



Hermite polynomials. Let us assume that the functions 
Ff n (9) are biorthogonal to G l k {9) in the interval [0, 6* max ], 



d0Gi(0)F4(0) = £ ro , fc , fc = 0,...,l. (45) 



The biorthogonality property ( fig) ensures that the inte- 
gral equation ( ^p| ) is fulfilled if the indices M; are con- 
structed in the same way as hi, Eq. (Q), where j p is 
replaced by m p + n p , 



N 



Mi = J^nip + ra p . 



(46) 



0=i 



Indeed, the integration over 6\ in (^) then selects cor- 
rect value of the sum mi+ni, subsequent integration over 
O2 fixes 7712 + n2 and so on. Notice also that the correct 
values of the differences nj — rrij are fixed by the exponen- 
tials exp[i(rij — mj)tpj] in (p9|). The sampling functions 
for multimode s-ordered moments are thus given by for- 
mula (39). 



The functions Fl n (9) biorthogonal to G k (8) have been 



discussed in 



One can construct them e.g. as linear 



combinations of G k {9), 



F l m {o) = Y, Al m k G l k {e). 



(47) 



k=0 



From the orthogonality conditions ( j45| ) one obtains a sys- 
tem of linear equations for the coefficients A l mk , which 
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can be solved for each m and /. If we choose # max = n, 
we can find simple analytical formulas for the functions 
F l 



k=0 



where (see Appendix for derivation) 



E 



m 
k-j 



(48) 



(-l) fe - J . (49) 



The functions F^(6) are plotted in Fig. 4. 

The sampling functions compensating imperfect detec- 
tion rj < 1, can be obtained from ( p9[ ) by making use of 
the simple replacement (Eoj), 



D mil (X,0,i/>;s,r]) = D r 



— ,6,ip;s 



v 



(50) 



This relation is particularly simple when normally or- 
dered moments are considered |16 . Inserting s = 1 into 
( |50| ) we have 

D ma (X, 0, -0; 1, 77) = r!- M / 2 D mn {X, 0, t/>; 1) . (51) 

Normally ordered moments do not contain any contribu- 
tion from vacuum fluctuations and they all vanish for a 
vacuum state. The experimental losses effectively reduce 
the value of normally ordered moment of M-th order by 
a factor r\ M l 2 . To compensate for imperfect detection, it 
suffices to use the ideal sampling function as if the detec- 
tion was perfect and then divide the result by rj M / 2 . 



B. Effect of aliasing and reconstruction limits 



In the experiment, the statistics w(X\ 8, if)) are mea- 
sured only at a certain finite number of angles 6j and 

(k) 

phases ip)j and the integration over dfl is replaced by a 
summation over finite number of discrete points (0,ip). 
This discretization imposes limits on the order of the re- 
constructed moments jL7j. 

Let us first consider the phases ipj. To simplify the 
discussion as much as possible, we restrict ourselves for a 
while to the single-mode case and sampling of symmet- 
rically ordered moments (s = 0, Weyl ordering). Let us 
further assume that the exact quadrature statistics are 
known for each of phases ij>^ k ' — kir/N^p. The sam- 
pling then reads, 



(at™a") sym = ^ 2 ( m +")/ 2 X(m,n) 

iV,' 



N,, 



k=l 



J^cxp ( i(n-m)— ) / dx 



:x m+n w I x, 



'4> / J — oo 



~n7, 



The formula (|52j) is a discrete Fourier transform in tjj. 
The Fourier series of the quadrature moments, 



dxx m+n w{x,iP) = 2-(" l+ ")/ 2 



m+n 

E 

fc=0 



m + n 



jm+n-fc fc\ i(m-j-n—2k)ip 



(53) 



contains either odd or even frequencies depending on the 
parity of m + n. The ./V^-point discrete Fourier transform 
(|52|) gives correct results only for sufficiently low mo- 
ments, because it cannot discriminate between exp[ifcf/>] 
and exp[i(k + 2N^)ip]. This phenomenon is called alias- 
ing |40f and it imposes an upper bound on the order of 
the reconstructed moment. When we substitute Fourier 
expansion (|5^) into Eq. (|52|), we find that m < and 
ri < Nij) must hold simultaneously. The same limitation 
obviously holds for any s-ordering and also for multimodc 
moment reconstruction with sampling functions D mll . In 
particular, 



rrij < , rij < N^ Jj , 



(54) 



(52) 



must be fulfilled, where N^. is the number of sampling 
points of the phase ipj . 

Let us proceed to the angles 6j . For a successful recon- 
struction, it is crucial to meet the biorthogonality condi- 
tions ( |45| ) where the integration is replaced by summation 
over Ng angles 9^ n \ 

Ng 

Y,F l m {e^)G l k {9^) = 5 mk , m,k = 0,...,l. (55) 

71=1 

If the condition (^5|) is violated due to discretization, then 
the reconstruction could be spoiled with large systematic 
error and the sampling would not yield reliable results. 
We shall prove below that the functions ( fl8l) fulfill the 
conditions ( pq ) provided that the sampling points are 
equidistant, IF™' = nn/Ng, and I < Ng holds. 

First of all we recall that the functions G l k (9), Eq. 
and F l m {6), Eq. (§§), can be expanded in finite Fourier 
series, with the highest component equal to I in both 
cases. Moreover, both functions contain only odd or only 
even Fourier components, depending on the parity of I. 
If Fj n (9) is expanded in Fourier series, then Eq. (|55|) be- 
comes a summation of several discrete Fourier transforms 
of G[(6) (we assume 9^ = kir/Ng). If I < Ng, then all 
discrete Fourier transforms yield the same results as the 
original integrations, and ([55|) holds exactly. The main 
advantage of the choice 6* max = it is now clear. It has al- 
lowed us to find analytical expressions for the functions 
FL which meet the discretized biorthogonality conditions 
©. We remark that the functions F l m n = F l m {9^) 
can also be constructed numerically by solving a system 
of Eqs. (pa) for a given set of sampling points 9^ [p0[. 
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Looking at formula ( |39| ) we find the limit on the order of 
reconstructed multimode moments, 



Mj < N e 



(56) 



We can conclude that if the quadrature statistics 
w(X] 8,ip) are measured with high accuracy, then sam- 
pling at finite number of points (6,ip) provides sufficient 
information for the successful reconstruction of certain s- 
ordered moments Cmn- If we use the sampling functions 
D mn and we want to reconstruct all moments of Mth 
order we have to measure w(X; 8, ip) at 



R(M, N) = (M + 1) 



2N-1 



(57) 



points (8, t/>) (we have the factor M + 1 for each of N 
phases ipj and N — 1 angles 9{). This number of sam- 
pling points is sufficient, but not necessary. The Mth 
order moments of TV-mode field can be parametrized by 
P{M, N) real numbers, where 



P(M, N) 



M 



■2N-1 
M 



(58) 



1 1 1 1 1 



It suffices to measure the statistics w(X;8,tp) at 
P(M,N) distinct points [8, if)). The appropriate sam- 
pling functions must be constructed numerically for a 
given set of sampling points (0,i/>) by inverting a sys- 
tem of linear equations which relates the moments C, 
to the moments of the quadrature statistics w(X;8,ip). 
This approach requires less sampling points because 
P(M,N) < R(M,N). Though many interesting ques- 
tions are related to this method, e.g. how to choose the 
points (9, ip), we do not deal with it in this paper in any 
more detail. 



Finally, we should note that the sampling functions are 
not unique. This is a general feature of optical homodyne 
tomography. An infinite number of functions exist whose 
average over w(x) 8, xjj) is zero for all physical quadrature 
distributions. These so-called null functions can be freely 
added to the above derived sampling functions. This free- 
dom of choice is exploited in adaptive homodyne tomog- 
raphy to find the sampling functions minimizing statisti- 
cal error for a given set of experimental data 13. 
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FIG. 5. Preparation of three-mode state. The squeezed 
vacuum in mode fei is mixed on beam splitters BSi and BS2 
with two vacua (modes 62 and 63), yielding the modes 01, 02, 
and as at the outputs. 



x 10" 



(C) 




FIG. 6. Reconstruction of the three-mode Q-function 
of a squeezed state prepared according to Fig. 5. A 
two-dimensional cut Q(a,a,a) through the six-dimensional 
phase space is plotted. Shown are surface (a) and contour (b) 
plots of the reconstructed quasidistribution and a difference 
AQ between reconstructed and exact Q- functions (c). 
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V. MONTE CARLO SIMULATIONS 

We have performed Monte Carlo simulations of multi- 
mode homodyne detection with a single LO and tested 
the performance of the sampling fnctions. Since the re- 
construction of multimode density matrix elements was 
already considered to relatively large extent in Ref. |32| , 
we focus here on the sampling of the multimode quasidis- 
tributions and s-ordered moments. The main purpose of 
this section is to illustrate the applicability of the above 
derived sampling functions and the feasibility of success- 
ful reconstruction of three-mode quantum state from an 
acceptably large amount of data. 

To be more specific, we consider three-mode squeezed 
state prepared according to Fig. 5. This state represents 
a simple but nontrivial example exhibiting nonclassical 
properties (squeezing). As depicted in Fig. 5, single- 
mode squeezed vacuum in mode 61 is mixed on two beam 
splitters BSi and BS2 with two vacua &2 and b%. The out- 
put modes 

ai = Tt 1 + 7i 62 ' 



(59) 



can then enter the multimode homodyne detector shown 
in Fig. 2. The transformation (|5{]) is unitary, thus pre- 
serving the canonical commutation relations. Moreover, 



b\ = bo cosh r + bt sinh r 



(60) 



where bo is annihilation operator of vacuum state and r is 
squeezing parameter. We assume r = 1 in the following. 

The reconstructed three-mode Q-hmctic-n j s shown in 
Fig. 6. In the computer simulation, we have sampled 



at 10 angles 9^ = /br/20, and phases njjf = 2vrfc/10, 
k = 1, 10. At each point (0,if>) the quadrature has 
been measured 50 times so that the total amount of ac- 
quired data is 5 x 10 6 . We have assumed a detection ef- 
ficiency 77 = 0.8 and we have used the loss-compensating 
sampling function (p0|). Q(ai, 0:2, 0:3) is a function in 
six-dimensional phase space, it is impossible to plot it as 
a whole and we must restrict ourselves to some lower- 
dimensional subspaces of the phase space. In Figure 6 
we show a two-dimensional cut Q(a, a, a). The recon- 
structed Q-function exhibits Gaussian shape characteris- 
tic for squeezed states. The squeezing is clearly reflected 
in the elliptic shape of the Q-function, as can be seen in 
the contour plot in Fig. 6(b). The reconstruction error 
can be judged from Fig. 6(c), which depicts the differ- 
ence AQ between exact and reconstructed Q-functions. 
The error is acceptably small and the reconstruction can 
be considered successful. 

Let us now proceed to sampling the multimode mo- 
ments. Again, we have assumed r\ = 0.8. We have 



,(fe) 



employed the loss-compensating sampling kernels (|5l| ) 
and the analytical functions F l m (9) given by Eq. (|48j). 
We have sampled at 10 different values of each angle 
6»| fe) = fc7r/10 and phase ^pf ] = kn/10, k = 1, . . . , 10. At 
each point (0,ip) 200 values of the quadrature X were 
recorded, which represents altogether 2 x 10 7 data. 

The reconstructed normally ordered moments (: n* :) 



and 



can be seen in Fig. 7. The gray bars display 



the exact values and allow for comparison with the sam- 
pled moments. The reconstructed moments are in very 
good agreement with the exact values. The sampling er- 
ror increases with the moment order and it is higher for 



(a!{) than for the factorial moments (:r ^ 
While (: n\ :) is still reconstructed with high accuracy, 
(af) is sampled with certain error. This can be explained 
by the necessity of sampling a high Fourier component 
exp(6i?/>i) in order to reconstruct (af). When we tried 
to sample moments of 10th or higher orders, the results 
suffered from very large systematic errors because we vi- 
olated the conditions ( |54| ) and (|56|). 

The moments (:n\ :) contain information on the pho- 
ton-number statistics of the mode a\ . From the sampled 
moments we can determine the Mandel Q-parameter for 
jth mode, 



Qj 



(:(An,) 2 :) (:nf.)-(n 3 ) 



( n j) 



(61) 



This parameter allows one to quickly distinguish between 
super-Poissonian (Qj > 0) and sub-Poissonian (Qj < 0) 
photon- number statistics. From the data shown in Fig. 
8 we have Qi w 1.25 and we find that the light in the 
mode cl\ exhibits super-Poissonian photon number statis- 
tics. The moments of the modes 02 and a.3 are the same 
as those of mode a\ because the squeezed vacuum b\ is 
equally split among the three modes a,j, c.f. Eq. (|5S|). 
The sampling works equally well for the modes 02 and 
as and the results are very similar to those displayed in 
Fig. 7. 





FIG. 7. Sampled moments of the mode ai. The empty solid 
bars show the reconstructed moments, the gray bars display 
exact values for comparison. Only real parts of the moments 
(0*) are shown. 
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FIG. 8. Sampled two-mode and three-mode moments. In 
Fig. (b), real parts of the complex moments are displayed. 

Having verified the feasibility of reconstruction of the 
single-mode moments we have finally sampled the mul- 
timode moments. Several results are shown in Fig. 8. 
Again, the low-order moments are reproduced with high 
accuracy, and the error increases with the moment order. 
It is worth noting that the photon number correlations 

(n* 1 n*f) a (62) 

can be sampled from phase averaged data. This is im- 
portant from the experimental point of view, because the 
sampling of moments (^32|) does not require stable rela- 
tive phase between the local oscillator and signal modes. 
All phases ipj can be fully randomized, e.g., by means of 
randomly driven piezoelectric modulators, and the homo- 
dyning then yields phase-averaged quadrature statistics 

In addition to the squeezed-vacuum state discussed 
here, we have also considered other quantum states, such 
as multimode coherent states and multimode squeezed 
coherent states. In all cases, the reconstruction pro- 
cedure worked well. The numerical simulations clearly 
demonstrate the feasibility of three-mode homodyne to- 
mography from w 10 7 recorded data. Of course, the 
number of necessary data inevitably increases with the 
number of modes. 



VI. CONCLUSIONS 

We have derived various important sampling functions 
for multimode homodyne tomography with a single lo- 
cal oscillator. Starting from the relation between mul- 
timode characteristic function and measured quadrature 



distribution we have found sampling functions for the s- 
parametrized quasidistributions with s < s v < 0. We 
have proved that the sampling function for Husimi qua- 
sidistribution is a generating function of the sampling 
functions / mn for density matrix elements in Fock basis 
p mn . The functions / mn were expressed as finite series 
of confluent hypergeometric functions. Finally, we have 
found the functions allowing for direct reconstruction of 
multimode s-ordered moments from the homodyne data. 
In all cases, loss-compensating sampling functions, appli- 
cable to a realistic experiment with detection efficiency 
77 < 1, have been provided. In order to test performance 
of the sampling functions we simulated homodyne de- 
tection of squeezed three-mode state and reconstructed 
the three-mode Q-function and several normally ordered 
moments. The reconstruction has shown very good re- 
sults for a detection efficiency i] = 0.8 and 10 7 sampled 
data, which is experimentally feasible. We emphasize 
that the multimode quantum state is reconstructed from 
the statistics of a class of single-mode quadratures. Only 
one homodyne detector is needed, which substantially 
simplifies the experiment. This method is particularly 
suitable for the measurement of ultrafast internal cor- 
relations of optical pulses or for the reconstruction of 
the quantum state of multimode single-frequency optical 
field. 
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APPENDIX: 

Here we derive the expression fl48[) f or the functions 
Ft^ifi). We insert the explicit form p2) of the function 
G(0) into @, multiply by a k {i(3) l - k and sum over k, 

Y,[ 0\aco S e) k (if3 S m9) l - k F} n (e)cie = a m (i(3) l - m 

(Al) 

The summation on the left-hand side is trivial and yields 




(acos9 + i(3sin9) l F^(e)de = a m (if3) l - m . (A2) 



In the next step we change variables, 6 = (a + (3)/2, 
7 = (a — (3)/2 and we have 
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{§e i» + je-^yp^dd = (7 + 6) m [i(6 - 7 )]' 



(A3) 



Now we set 6 = 1, differentiate (A3) /c-times with respect 
to 7 and then set 7 = 0. After little algebra we arrive at 



,i(l-2k)6 P ! 



FL(0)d9 



{l-k)l d k 



[(l+7) m (l-7) 



I— ml 



(A4) 



7=0 



Now we assume that the function F l m can be written in 
terms of finite Fourier series, 



(A5) 



and insert this expansion into (A4). After integration on 
the left -han d side and differentiation on the right-hand 
side of (A4) we find 



EL 



' "' // \ 1 / rn\ (I — m 



E 

j=0 



3 J \n-J 



: (A6) 



and we have derived the formulas (Eq) and (p 
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